Comparison between trajectories in the Gmnc WT/KO

Load libraries

library(Seurat)
library(princurve)
library(Revelio)
library(monocle)
library(gprofiler2)
library(seriation)
library(Matrix)
library(dplyr)
library(RColorBrewer)
library(fungible)
library(ggplot2)
library(ggExtra)
library(cowplot)
library(wesanderson)

#Set ggplot theme as classic
theme_set(theme_classic())

Pseudotime in the WT

WT <- readRDS("../QC.filtered.clustered.cells.RDS")

Extract CR and CP neurons

WT$Lineage <- sapply(WT$Cell_ident,
                              FUN = function(x) {
                                if (x %in% c("Cajal-Retzius_neurons", "Hem")) {
                                  x = "Cajal-Retzius_neurons"
                                } else if (x %in% c("Pallial_neurons", "Medial_pallium")) {
                                  x = "Pallial_neurons"
                                } else if(x %in% c("ChP", "ChP_progenitors")) {
                                   x =  "Choroid_Plexus"
                                } else {
                                  x = "other"
                                  }
                              })

Neurons.data <-  subset(WT, idents = c("Cajal-Retzius_neurons", "Pallial_neurons"))

DimPlot(Neurons.data,
        reduction = "spring",
        pt.size = 1,
        cols =  c("#cc391b","#026c9a")
        ) + NoAxes()

rm(WT)

Fit principale curve on the two lineages

Cajal-Retzius cells

Trajectories.Hem <- Neurons.data@meta.data %>%
                    select("Barcodes", "nUMI", "Spring_1", "Spring_2", "AP_signature1","BP_signature1", "EN_signature1", "LN_signature1", "Lineage") %>%
                    filter(Lineage == "Cajal-Retzius_neurons")

fit <- principal_curve(as.matrix(Trajectories.Hem[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = .7,
                       stretch=0)

## Starting curve---distance^2: 45804778678
## Iteration 1---distance^2: 27732113
## Iteration 2---distance^2: 27728318

#The principal curve smoothed
Hem.pc.line <- as.data.frame(fit$s[order(fit$lambda),]) 

#Pseudotime score
Trajectories.Hem$PseudotimeScore <- fit$lambda/max(fit$lambda)

if (cor(Trajectories.Hem$PseudotimeScore, Neurons.data@assays$SCT@data['Hmga2', Trajectories.Hem$Barcodes]) > 0) {
  Trajectories.Hem$PseudotimeScore <- -(Trajectories.Hem$PseudotimeScore - max(Trajectories.Hem$PseudotimeScore))
}

Pallial neurons

Trajectories.Pallial <- Neurons.data@meta.data %>%
                        select("Barcodes", "nUMI", "Spring_1", "Spring_2", "AP_signature1","BP_signature1", "EN_signature1", "LN_signature1", "Lineage") %>%
                        filter(Lineage == "Pallial_neurons")

fit <- principal_curve(as.matrix(Trajectories.Pallial[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = .7,
                       stretch=0)

## Starting curve---distance^2: 26984853690
## Iteration 1---distance^2: 22153700
## Iteration 2---distance^2: 22179462
## Iteration 3---distance^2: 22180297

#The principal curve smoothed
Pallial.pc.line <- as.data.frame(fit$s[order(fit$lambda),])

#Pseudotime score
Trajectories.Pallial$PseudotimeScore <- fit$lambda/max(fit$lambda)

if (cor(Trajectories.Pallial$PseudotimeScore, Neurons.data@assays$SCT@data['Hmga2', Trajectories.Pallial$Barcodes]) > 0) {
  Trajectories.Pallial$PseudotimeScore <- -(Trajectories.Pallial$PseudotimeScore - max(Trajectories.Pallial$PseudotimeScore))
}

Combine the two trajectories’ data

Trajectories.neurons.WT <- rbind(Trajectories.Pallial, Trajectories.Hem)

cols <- brewer.pal(n =11, name = "Spectral")

ggplot(Trajectories.neurons.WT, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=PseudotimeScore), size=2, shape=16) + 
  scale_color_gradientn(colours=rev(cols), name='Speudotime score') +
  geom_line(data=Pallial.pc.line, color="#026c9a", size=0.77) +
  geom_line(data=Hem.pc.line, color="#cc391b", size=0.77)

Plot pan-neuronal genes along this axis

Neurons.data <- NormalizeData(Neurons.data, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

# Neurog2
p1 <- FeaturePlot(object = Neurons.data,
            features = c("Neurog2"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

Trajectories.neurons.WT$Neurog2 <- Neurons.data@assays$RNA@data["Neurog2", Trajectories.neurons.WT$Barcodes]

p2 <- ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore, y= Neurog2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Tbr1 
p3 <- FeaturePlot(object = Neurons.data ,
            features = c("Tbr1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()
Trajectories.neurons.WT$Tbr1 <- Neurons.data@assays$RNA@data["Tbr1", Trajectories.neurons.WT$Barcodes]

p4 <- ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore, y= Tbr1)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Mapt 
p5 <- FeaturePlot(object = Neurons.data ,
            features = c("Mapt"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

Trajectories.neurons.WT$Mapt <- Neurons.data@assays$RNA@data["Mapt", Trajectories.neurons.WT$Barcodes]

p6 <- ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore, y= Mapt)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)

Shift Pseudotime in both lineage

Since we observe the first 25% of both trajectories are occupied by few, likely progenitor cells, we shift this cell along the axis

Pseudotime.intervals <- Trajectories.neurons.WT%>%
                          select(Lineage, PseudotimeScore) %>%
                          mutate(Pseudotime.bins = cut(Trajectories.neurons.WT$PseudotimeScore, seq(0, max(Trajectories.neurons.WT$PseudotimeScore) + 0.05, 0.05), dig.lab = 2, right = FALSE)) %>%
                          group_by(Lineage, Pseudotime.bins) %>%
                          summarise(n=n())

ggplot(Pseudotime.intervals, aes(x=Pseudotime.bins, y=n, fill=Lineage)) +
        geom_bar(stat = "identity", width = 0.90) +
        theme(axis.text.x = element_text(angle = 45, hjust=1))+
        scale_fill_manual(values= c("#cc391b", "#026c9a"))

score <- sapply(Trajectories.neurons.WT$PseudotimeScore,
                FUN = function(x) if (x <= 0.2) {x= 0.2} else { x=x })

Trajectories.neurons.WT$PseudotimeScore.shifted <- (score - min(score)) / (max(score) - min(score))

# Neurog2
p1 <- FeaturePlot(object = Neurons.data ,
            features = c("Neurog2"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore.shifted, y= Neurog2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Tbr1 
p3 <- FeaturePlot(object = Neurons.data ,
            features = c("Tbr1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p4 <- ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore.shifted, y= Tbr1)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Mapt 
p5 <- FeaturePlot(object = Neurons.data ,
            features = c("Mapt"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p6 <- ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore.shifted, y= Mapt)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)

ggplot(Trajectories.neurons.WT, aes(x= PseudotimeScore.shifted, y= nUMI/10000)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

Hem.data <- readRDS("../QC.filtered.clustered.cells.RDS")

Trajectories.neurons.WT$Pseudotime <- Trajectories.neurons.WT$PseudotimeScore.shifted

Neuro.trajectories_WT <- CreateSeuratObject(counts = Hem.data@assays$RNA@data[, Trajectories.neurons.WT$Barcodes],
                                            meta.data = Trajectories.neurons.WT)

spring <- as.matrix(Neuro.trajectories_WT@meta.data %>% select("Spring_1", "Spring_2"))
  
Neuro.trajectories_WT[["spring"]] <- CreateDimReducObject(embeddings = spring, key = "Spring_", assay = DefaultAssay(Neuro.trajectories_WT))

p1 <- FeaturePlot(object = Neuro.trajectories_WT,
            features = "Pseudotime",
            pt.size = 0.5,
            cols = rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- DimPlot(object = Neuro.trajectories_WT,
        group.by = "Lineage",
        pt.size = 0.5,
        reduction = "spring",
        cols =  c("#cc391b", "#026c9a")) & NoAxes()

p1 + p2

rm(list = ls()[!ls() %in% c("Neuro.trajectories_WT")])

Pseudotime in the KO

KO <- readRDS("./GmncKO.cells.RDS") %>% subset(idents = c(6:9), invert = T)

Compute differentiation states scores

AP

APgenes <- c("Rgcc", "Sparc", "Hes5","Hes1", "Slc1a3",
             "Ddah1", "Ldha", "Hmga2","Sfrp1", "Id4",
             "Creb5", "Ptn", "Lpar1", "Rcn1","Zfp36l1",
             "Sox9", "Sox2", "Nr2e1", "Ttyh1", "Trip6")

KO <- AddModuleScore(KO,
                     features = list(APgenes),
                     name = "AP_signature")

BP

BPgenes <- c("Eomes", "Igsf8", "Insm1", "Elavl2", "Elavl4",
             "Hes6","Gadd45g", "Neurog2", "Btg2", "Neurog1")

KO <- AddModuleScore(KO,
                     features = list(BPgenes),
                     name = "BP_signature")

EN

ENgenes <- c("Mfap4", "Nhlh2", "Nhlh1", "Ppp1r14a", "Nav1",
             "Neurod1", "Sorl1", "Svip", "Cxcl12", "Tenm4",
             "Dll3", "Rgmb", "Cntn2", "Vat1")

KO <- AddModuleScore(KO,
                     features = list(ENgenes),
                     name = "EN_signature")

LN

LNgenes <- c("Snhg11", "Pcsk1n", "Mapt", "Ina", "Stmn4",
             "Gap43", "Tubb2a", "Ly6h","Ptprd", "Mef2c")

KO <- AddModuleScore(KO,
                     features = list(LNgenes),
                     name = "LN_signature")

FeaturePlot(object = KO,
            features = c("AP_signature1", "BP_signature1",
                              "EN_signature1", "LN_signature1"),
            pt.size = 0.75,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "spring",
            order = T) & NoAxes() & NoLegend()

Group cells in Pallial or CR lineage

KO$Lineage <- sapply(KO$Cell.ident,
                              FUN = function(x) {
                                if (x %in% c("Neuron_prob.2", "Hem")) {
                                  x = "Cajal-Retzius_neurons"
                                } else if (x %in% c("Neuron_prob.3", "Medial_pallium")) {
                                  x = "Pallial_neurons"
                                } else {
                                  x = "other"
                                  }
                              })

DimPlot(KO,
        reduction = "spring",
        group.by = "Lineage",
        pt.size = 0.5,
        cols =  c("#cc391b","#969696","#026c9a")
        ) + NoAxes()

Fit principale curve on the two lineages

Neurons.data <-  subset(KO,  subset = Lineage %in% c("Cajal-Retzius_neurons", "Pallial_neurons") & Cell.ident %in% c("Neuron_prob.2", "Neuron_prob.3"))

DimPlot(Neurons.data ,
        reduction = "spring",
        group.by = "Lineage",
        pt.size = 1,
        cols =  c("#cc391b","#026c9a")
        ) + NoAxes()

fit <- principal_curve(as.matrix(Neurons.data@meta.data[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 1,
                       stretch=0)

## Starting curve---distance^2: 76242448164
## Iteration 1---distance^2: 48308567
## Iteration 2---distance^2: 49237931
## Iteration 3---distance^2: 50119026
## Iteration 4---distance^2: 51105768
## Iteration 5---distance^2: 51819999
## Iteration 6---distance^2: 52369574
## Iteration 7---distance^2: 52732675
## Iteration 8---distance^2: 52936995
## Iteration 9---distance^2: 53060989
## Iteration 10---distance^2: 53117718

#Pseudotime score
PseudotimeScore <- fit$lambda/max(fit$lambda)

if (cor(PseudotimeScore, Neurons.data@assays$SCT@data['Hmga2', ]) > 0) {
  Neurons.data$PseudotimeScore <- -(PseudotimeScore - max(PseudotimeScore))
}

cols <- brewer.pal(n =11, name = "Spectral")

ggplot(Neurons.data@meta.data, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=PseudotimeScore), size=2, shape=16) + 
  scale_color_gradientn(colours=rev(cols), name='Pseudotime score')

Plot pan-neuronal genes along this axis

Neurons.data <- NormalizeData(Neurons.data, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

Trajectories.neurons <- Neurons.data@meta.data %>% select(Barcodes, Spring_1, Spring_2,
                                                          AP_signature1, BP_signature1, EN_signature1, LN_signature1,
                                                          Lineage, PseudotimeScore)

# Neurog2
p1 <- FeaturePlot(object = Neurons.data,
            features = c("Neurog2"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

Trajectories.neurons$Neurog2 <- Neurons.data@assays$RNA@data["Neurog2", Trajectories.neurons$Barcodes]

p2 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Neurog2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Tbr1 
p3 <- FeaturePlot(object = Neurons.data ,
            features = c("Tbr1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()
Trajectories.neurons$Tbr1 <- Neurons.data@assays$RNA@data["Tbr1", Trajectories.neurons$Barcodes]

p4 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Tbr1)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Mapt 
p5 <- FeaturePlot(object = Neurons.data ,
            features = c("Mapt"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

Trajectories.neurons$Mapt <- Neurons.data@assays$RNA@data["Mapt", Trajectories.neurons$Barcodes]

p6 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Mapt)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)

Shift pseudotime score

score <- sapply(Trajectories.neurons$PseudotimeScore,
                FUN = function(x) if (x <= 0.15) {x= 0.15} else { x=x })

Trajectories.neurons$PseudotimeScore.shifted <- (score - min(score)) / (max(score) - min(score))

# Neurog2
p1 <- FeaturePlot(object = Neurons.data ,
            features = c("Neurog2"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Neurog2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Tbr1 
p3 <- FeaturePlot(object = Neurons.data ,
            features = c("Tbr1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p4 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Tbr1)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Mapt 
p5 <- FeaturePlot(object = Neurons.data ,
            features = c("Mapt"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p6 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Mapt)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)

Trajectories.neurons$nUMI <- Neurons.data$nCount_RNA

ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= nUMI/10000)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

Subset the full Seurat object

KO <- readRDS("./GmncKO.cells.RDS") %>% subset(idents = c(6:9), invert = T)

Trajectories.neurons$Pseudotime <- Trajectories.neurons$PseudotimeScore.shifted

Neuro.trajectories_KO <- CreateSeuratObject(counts = KO@assays$RNA@data[, Trajectories.neurons$Barcodes],
                                            meta.data = Trajectories.neurons)

spring <- as.matrix(Neuro.trajectories_KO@meta.data %>% select("Spring_1", "Spring_2"))
  
Neuro.trajectories_KO[["spring"]] <- CreateDimReducObject(embeddings = spring, key = "Spring_", assay = DefaultAssay(Neuro.trajectories_KO))

p1 <- FeaturePlot(object = Neuro.trajectories_KO,
            features = "Pseudotime",
            pt.size = 0.5,
            cols = rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- DimPlot(object = Neuro.trajectories_KO,
        group.by = "Lineage",
        pt.size = 0.5,
        reduction = "spring",
        cols =  c("#cc391b", "#026c9a")) & NoAxes()

p1 + p2

Cell euclidian distance along pseudotime

rm(list = ls()[!ls() %in% c("Neuro.trajectories_WT","Neuro.trajectories_KO")])
gc()

##            used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells  3515984 187.8    6012948  321.2   6012948  321.2
## Vcells 26894381 205.2  753491556 5748.7 890495468 6794.0

Group cells in 10 bins along pseudotime and average expression by lineage and genotype

TFs <- read.table("TF.csv", sep = ";")[,1]

WT_KO <- merge(x = Neuro.trajectories_WT, y = Neuro.trajectories_KO)

WT_KO <- NormalizeData(WT_KO, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

WT_KO <-FindVariableFeatures(WT_KO, selection.method = "disp", nfeatures = 3000, assay = "RNA")

WT_KO  <- ScaleData(WT_KO, vars.to.regress = "nCount_RNA", features = TFs)

WT_KO$Genotype <- sapply(WT_KO$orig.ident,
                              FUN = function(x) {
                                if (x %in% c("Hem1", "Hem2")) {
                                  x = "WT"
                                } else {
                                  x = "KO"
                                  }
                              })
nbreaks <- 5
WT_KO$Pseudotime.bins <- cut(WT_KO$Pseudotime, breaks = nbreaks , labels = 1:nbreaks)

# WT Cajal-Retzius
WT_CR <- subset(WT_KO, subset = Genotype == "WT" & Lineage == "Cajal-Retzius_neurons")
  
WT_CR.averages <- AverageExpression(WT_CR,
                                    features = TFs,
                                    group.by = "Pseudotime.bins",
                                    slot = "scale.data")

# KO Cajal-Retzius
KO_CR <- subset(WT_KO, subset = Genotype == "KO" & Lineage == "Cajal-Retzius_neurons")

KO_CR.averages <- AverageExpression(KO_CR,
                                    features = TFs,
                                    group.by = "Pseudotime.bins",
                                    slot = "scale.data")
  
# WT Pallial neurons
WT_Pal <- subset(WT_KO, subset = Genotype == "WT" & Lineage == "Pallial_neurons")

WT_Pal.averages <- AverageExpression(WT_Pal,
                                    features = TFs,
                                    group.by = "Pseudotime.bins",
                                    slot = "scale.data")
  
# KO Pallial neurons  
KO_Pal <- subset(WT_KO, subset = Genotype == "KO" & Lineage == "Pallial_neurons")

KO_Pal.averages <- AverageExpression(KO_Pal,
                                    features = TFs,
                                    group.by = "Pseudotime.bins",
                                    slot = "scale.data")

Compute correlation across conditions

ggplot(WT_CR@meta.data, aes(Spring_1, Spring_2, color= Pseudotime.bins)) + geom_point()

### Distance to KO pallial neurons

Smoothed.point.cor <- pdist::pdist(X= t(WT_CR.averages$RNA),
                                   Y= t(KO_Pal.averages$RNA))

WT_CR.dist <- diag(as.matrix(Smoothed.point.cor))

#
Smoothed.point.cor <- pdist::pdist(X= t(KO_CR.averages$RNA),
                                   Y= t(KO_Pal.averages$RNA))


KO_CR.dist <- diag(as.matrix(Smoothed.point.cor))

#
Smoothed.point.cor <- pdist::pdist(X= t(WT_Pal.averages$RNA),
                                   Y= t(KO_Pal.averages$RNA))

WT_Pal.dist <- diag(as.matrix(Smoothed.point.cor))

data <- data.frame(Distance = c(WT_CR.dist, KO_CR.dist, WT_Pal.dist),
                   Trajectory = rep(c("WT_CR", "KO_CR", "WT_Pal"), each=nbreaks),
                   Pseudotime.bin = factor(rep(paste0("Bin_",1:nbreaks), 3),levels = paste0("Bin_",1:nbreaks)))

ggplot(data, aes(x= Pseudotime.bin, y= Distance, color= Trajectory, group = Trajectory)) +
  geom_line() +
  geom_point() + ggtitle("Distance to KO pallial neurons")

### Distance to WT pallial neurons

Smoothed.point.cor <- pdist::pdist(X= t(WT_CR.averages$RNA),
                                   Y= t(WT_Pal.averages$RNA))

WT_CR.dist <- diag(as.matrix(Smoothed.point.cor))

#
Smoothed.point.cor <- pdist::pdist(X= t(KO_CR.averages$RNA),
                                   Y= t(WT_Pal.averages$RNA))


KO_CR.dist <- diag(as.matrix(Smoothed.point.cor))

#
Smoothed.point.cor <- pdist::pdist(X= t(KO_Pal.averages$RNA),
                                   Y= t(WT_Pal.averages$RNA))

KO_Pal.dist <- diag(as.matrix(Smoothed.point.cor))

data <- data.frame(Distance = c(WT_CR.dist, KO_CR.dist, KO_Pal.dist),
                   Trajectory = rep(c("WT_CR", "KO_CR", "KO_Pal"), each=nbreaks),
                   Pseudotime.bin = factor(rep(paste0("Bin_",1:nbreaks), 3),levels = paste0("Bin_",1:nbreaks)))

ggplot(data, aes(x= Pseudotime.bin, y= Distance, color= Trajectory, group = Trajectory)) +
  geom_line() +
  geom_point() + ggtitle("Distance to WT pallial neurons")

Cosine to KO pallial neurons

Smoothed.point.cor <- cosMat(WT_CR.averages$RNA,
                             KO_Pal.averages$RNA)

WT_CR.dist <- diag(Smoothed.point.cor$cosine)

#
Smoothed.point.cor <- cosMat(KO_CR.averages$RNA,
                             KO_Pal.averages$RNA)


KO_CR.dist <- diag(Smoothed.point.cor$cosine)

#
Smoothed.point.cor <- cosMat(WT_Pal.averages$RNA,
                             KO_Pal.averages$RNA)

KO_Pal.dist <- diag(Smoothed.point.cor$cosine)

data <- data.frame(Distance = c(WT_CR.dist, KO_CR.dist, KO_Pal.dist),
                   Trajectory = rep(c("WT_CR", "KO_CR", "WT_Pal"), each=nbreaks),
                   Pseudotime.bin = factor(rep(paste0("Bin_",1:nbreaks), 3),levels = paste0("Bin_",1:nbreaks)))

ggplot(data, aes(x= Pseudotime.bin, y= Distance, color= Trajectory, group = Trajectory)) +
  geom_line() +
  geom_point() + ggtitle("Distance to KO pallial neurons")

Session Info

#date
format(Sys.time(), "%d %B, %Y, %H,%M")

## [1] "03 juin, 2022, 18,42"

#Packages used
sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] splines   stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] wesanderson_0.3.6   cowplot_1.1.1       ggExtra_0.9        
##  [4] fungible_2.0        RColorBrewer_1.1-2  dplyr_1.0.7        
##  [7] seriation_1.3.1     gprofiler2_0.2.1    monocle_2.22.0     
## [10] DDRTree_0.1.5       irlba_2.3.3         VGAM_1.1-5         
## [13] ggplot2_3.3.5       Biobase_2.54.0      BiocGenerics_0.40.0
## [16] Matrix_1.4-1        Revelio_0.1.0       princurve_2.1.6    
## [19] SeuratObject_4.0.4  Seurat_4.0.5       
## 
## loaded via a namespace (and not attached):
##   [1] plyr_1.8.6            igraph_1.2.11         GPArotation_2014.11-1
##   [4] lazyeval_0.2.2        densityClust_0.3      listenv_0.8.0        
##   [7] scattermore_0.7       fastICA_1.2-3         digest_0.6.29        
##  [10] foreach_1.5.1         htmltools_0.5.2       viridis_0.6.2        
##  [13] fansi_0.5.0           magrittr_2.0.2        tensor_1.5           
##  [16] cluster_2.1.3         ROCR_1.0-11           limma_3.50.0         
##  [19] globals_0.14.0        matrixStats_0.61.0    docopt_0.7.1         
##  [22] spatstat.sparse_2.0-0 pdist_1.2.1           colorspace_2.0-2     
##  [25] ggrepel_0.9.1         xfun_0.28             sparsesvd_0.2        
##  [28] crayon_1.4.2          jsonlite_1.7.2        spatstat.data_2.1-0  
##  [31] survival_3.2-13       zoo_1.8-9             iterators_1.0.13     
##  [34] glue_1.5.1            polyclip_1.10-0       registry_0.5-1       
##  [37] gtable_0.3.0          leiden_0.3.9          future.apply_1.8.1   
##  [40] abind_1.4-5           scales_1.1.1          mvtnorm_1.1-3        
##  [43] pheatmap_1.0.12       DBI_1.1.1             miniUI_0.1.1.1       
##  [46] Rcpp_1.0.8            viridisLite_0.4.0     xtable_1.8-4         
##  [49] clue_0.3-60           reticulate_1.22       spatstat.core_2.3-1  
##  [52] htmlwidgets_1.5.4     httr_1.4.2            FNN_1.1.3            
##  [55] ellipsis_0.3.2        ica_1.0-2             farver_2.1.0         
##  [58] pkgconfig_2.0.3       sass_0.4.0            uwot_0.1.10          
##  [61] deldir_1.0-6          utf8_1.2.2            labeling_0.4.2       
##  [64] tidyselect_1.1.1      rlang_0.4.12          reshape2_1.4.4       
##  [67] later_1.3.0           munsell_0.5.0         tools_4.2.0          
##  [70] generics_0.1.1        ggridges_0.5.3        evaluate_0.14        
##  [73] stringr_1.4.0         fastmap_1.1.0         yaml_2.2.1           
##  [76] goftest_1.2-3         knitr_1.36            fitdistrplus_1.1-6   
##  [79] purrr_0.3.4           RANN_2.6.1            pbapply_1.5-0        
##  [82] future_1.23.0         nlme_3.1-153          mime_0.12            
##  [85] slam_0.1-49           compiler_4.2.0        plotly_4.10.0        
##  [88] png_0.1-7             spatstat.utils_2.2-0  tibble_3.1.6         
##  [91] bslib_0.3.1           stringi_1.7.6         highr_0.9            
##  [94] lattice_0.20-45       HSMMSingleCell_1.14.0 vctrs_0.3.8          
##  [97] pillar_1.6.4          lifecycle_1.0.1       spatstat.geom_2.3-0  
## [100] combinat_0.0-8        lmtest_0.9-39         jquerylib_0.1.4      
## [103] RcppAnnoy_0.0.19      data.table_1.14.2     httpuv_1.6.3         
## [106] patchwork_1.1.1       R6_2.5.1              promises_1.2.0.1     
## [109] TSP_1.1-11            KernSmooth_2.23-20    gridExtra_2.3        
## [112] nleqslv_3.3.2         parallelly_1.29.0     codetools_0.2-18     
## [115] MASS_7.3-57           assertthat_0.2.1      withr_2.4.3          
## [118] qlcMatrix_0.9.7       sctransform_0.3.2     mgcv_1.8-40          
## [121] parallel_4.2.0        grid_4.2.0            rpart_4.1.16         
## [124] Rcsdp_0.1.57.2        tidyr_1.1.4           rmarkdown_2.11       
## [127] Rtsne_0.15            shiny_1.7.1

---

1. Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩︎